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FILTERED BACKPROJECTION ALGORITHMS FOR COMPTON 
CAMERAS IN NUCLEAR MEDICINE 

[01] This application claims priority to provisional U.S. Application. No. 60/458,131, filed 
March 27, 2003. 

FIELD OF THE INVENTION 

[02] The present invention relates to nuclear medicine for reconstructing a source distribution 
from Compton camera data in the application of nuclear medicine. In particular, the 
source distribution corresponds to a radiopharmaceutical distribution within a patient. 

BACKGROUND OF THE INVENTION 

[03] In nuclear medicine, a radiopharmaceutical is injected into a patient and accumulates in 
tissues of physiological interest (e.g., cancer tumors or heart muscle). The gamma-ray 
emissions from this radiopharmaceutical are then observed and the local density of 
radiopharmaceutical within the patient reconstructed. Two clinical imaging techniques, 
single photon emission computed tomography (SPECT) and positron emission 
tomography (PET), are based on this simple idea and are currently available 
commercially. Both these imaging techniques use filtered backprojection (FBP) 
algorithms for the reconstruction of the three-dimensional (3D) distribution of 
radiopharmaceuticals within the patient. The advantages and limitations of SPECT and 
PET have been reviewed extensively in the literature and will not be further discussed 
here. It is generally acknowledged, however, that one of the maj or limitations of SPECT 
is the loss of counts due to the collimator in the gamma-ray camera. Typically only one 
gamma ray is observed for every ten thousand emissions, the vast majority of these 
gamma rays are absorbed in the collimator. In 1974 Todd, Nightingale, and Everett 
proposed a new type of gamma camera that would replace the collimator with a method 
of "electronic collimation" based on Compton scattering. Compton kinematics implies a 
relationship between the angle of scatter and the energy deposited in a detector. 
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Knowledge of this energy and, therefore, the angle of scatter can be used to estimate the 
direction of an incident gamma-ray. The success of this strategy is notably demonstrated 
by the gamma-ray telescopes currently surveying the sky for sources of gamma emissions. 
In medical applications, however, the patient is located at a finite distance, so that the 
interpretation of the angular information provided by Compton events is complicated by 
huge parallax effects. In landmark research during the 1980s and 1990s, Manbir Singh 
and his group constructed and tested a prototype Compton camera. That research 
revealed both the promise and difficulties associated with Compton cameras. Many of 
the technical difficulties associated with Compton cameras (e.g., limited energy 
resolution, restrictive detector geometries, and events involving multiple scatters) are 
currently being remedied by advances in detector technology. Among the most 
significant of these problems is the 3D reconstruction of the radiopharmaceutical 
distribution from the data provided by a Compton camera. 

The reconstruction of 3D images from Compton camera data poses a formidable 
mathematical problem. Unlike collimated cameras that reveal the direction of each 
gamma-ray, the incident direction of each event in a Compton camera is restricted to a 
cone. This geometric difference completely changes the nature of the reconstruction 
problem. In his original studies, Singh in collaboration with Hebert and Leahy 
implemented a maximum likelihood (ML) reconstruction (and tested other iterative 
methods). In recent years numerous alternative algorithms have been proposed and 
studied. Roughly speaking, these algorithms can be divided into iterative and analytic 
inversion methods. Both types of algorithm require the construction of a mathematical 
model that describes the Compton camera data in terms of the source distribution of 
radiopharmaceutical substance. Generally, the iterative methods can accommodate more 
intricate and accurate models, whereas, the analytic inversions require simpler models 
that are amenable to direct mathematical analysis. Most experimental groups working on 
Compton cameras have used iterative algorithms similar to the ML techniques of Singh et 
al. However, the reason for preferring iterative algorithms was not the accuracy of the 
mathematical model. After all, FBP algorithms (that are prototypes of analytic inversion) 
are generally used in clinical SPECT and PET despite the existence of ML algorithms for 
more accurate models. Basically, the iterative algorithms dominate Compton 
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reconstruction because, until the late 1990s, no analytic algorithm existed for Compton 
cameras. 

An algorithm for Compton cameras developed by Cree and Bones restricts the acceptable 
data as much as would a standard collimator — and is, therefore, unsuitable for practical 
applications. Shortly thereafter, Basko, Gullberg and Zeng produced and patented an 
algorithm (for the same mathematical model) that removed all restrictions on the 
Compton data. This algorithm explicitly rebins the Compton data into a form that leads 
to a 2D Radon transform (and, therefore, FBP). This rebinning procedure, which is the 
major drawback of the technique, overcomes the basic problem of Compton camera 
reconstruction; namely, the cone of incident gamma-ray directions. The rebinning 
implements a complex summation over spherical harmonics that arises from the cone of 
acceptance. Other researchers developed analytic-inversion algorithms that explicitly 
replace integrals over the acceptance cone with a summation over spherical harmonics. 
For example, Parra and Tomitani and Hirasawa retain the spherical harmonic expansions 
— truncating the series to a finite number of terms for numerical implementation. 
Although Parra and Tomitani and Hirasawa use a slightly different mathematical model 
than Basko et al, they face the same problem; namely, the additional degree of freedom 
introduced by the acceptance cone of the Compton camera. In reading the papers of Parra 
and Tomitani and Hirasawa, one is reminded of the initial papers of Cormack in which 
the Radon transform was rediscovered in terms of circular harmonics. One recalls that, 
when the problem was properly reformulated, the circular harmonics disappeared and the 
corresponding FBP algorithms emerged. 

While these analytic algorithms were being developed, other groups examined alternative 
strategies. In the process of designing and building the C-SPRINT system, the group at 
the University of Michigan engaged in a broad-based campaign to establish Compton 
cameras as a viable imaging technology. Not only were simulations and sophisticated 
statistical analysis used for camera design, but the reconstruction problem was attacked 
using a variety of techniques. As a leading research institution in ML and statistical 
reconstruction, the University of Michigan implemented the most sophisticated iterative 
algorithms currently available. Furthermore, non-iterative reconstruction algorithms were 
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also examined. Sauve et al studied matrix inversion techniques for a model of their 
system. The geometrical symmetries of the system were exploited in an attempt to handle 
the otherwise intractably large set of equations. In another approach, that is relevant to 
the research reported in this paper, Wilderman et al developed the software for list-mode 
back-projection of Compton scatter events into the appropriate cone geometries. In this 
research, they were joined by Rohe et al who also developed backprojection algorithms 
for the cone geometry of Compton cameras. Unfortunately, simple backprojection onto 
the appropriate cones (without filtering) does not produce accurate reconstruction of the 
source distribution. 

With the current art, reconstruction of a radiopharmaceutical source distribution typically 
utilizes iterative algorithms. The current art also provides analytic-inversion techniques 
that may be inherently faster to execute but are typically deficient with respect to relative 
immunity to statistical noise and the desired accuracy of reconstructing the source 
distribution. Hence, there is a real need to develop methods and apparatus that utilize 
analytic-inversion techniques while providing accuracy in reconstructing a 
radiopharmaceutical source distribution within a patient in order to properly diagnose the 
patient. 

BRIEF SUMMARY OF THE INVENTION 

The present invention provides methods and apparatuses for acquisition and processing of 
the observed data from a Compton camera in order to construct a three-dimensional 
image of a radiopharmaceutical source distribution within a patient. No intermediate two- 
dimensional images are formed. 

Another aspect of the invention provides apparatuses that include a Compton camera for 
obtaining observed data from a radiopharmaceutical source distribution. The observed 
data is analyzed by a processor in order to construct a three-dimensional image 
representing the source distribution. The three-dimensional image may be displayed for a 
clinician. Moreover, the clinician may provide input data for configuring the processor for 
processing the observed data. 
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[10] With another aspect of the invention, an idealized mathematical model expresses the 
observed Compton camera data in terms of an integral over the source distribution. An 
exact analytic inversion is then found for this idealized model. The new analytic 
solutions arise from a generalization of the integral equation used to model the Compton 
camera. The kernel of the integral equation is modified by the introduction of an index p 
that describes the effect of source distance from the camera. 

BRIEF DESCRIPTION OF THE DRAWINGS 

[11] Figure 1 shows a basic Compton cones geometry in accordance with an embodiment of 
the invention. 

[12] Figure 2A shows a rectangular detector in accordance with an embodiment of the 
invention. 

[13] Figure 2B shows a square prism detector in accordance with an embodiment of the 
invention. 

[14] Figure 2C shows a cylindrical geometry in accordance with an embodiment of the 
invention. 

[15] Figure 2D shows a circular geometry in accordance with an embodiment of the invention. 

[16] Figure 2E shows a spherical detector in accordance with an embodiment of the invention. 

[17] Figure 3 shows alternative detector D2 geometries in accordance with an embodiment of 
the invention. 

[18] Figure 4 A shows alternative scanning geometries for a linear scan in accordance with an 
embodiment of the invention. 

[19] Figure 4B shows alternative scanning geometries for an orbital scan in accordance with 
an embodiment of the invention. 

[20] Figure 5 shows a conceptual flowchart for filtered backprojection (FBP) algorithms in 
accordance with an embodiment of the invention. 

[21] Figure 6 shows a detailed flowchart of filtered backprojection (FBP) algorithms in 
accordance with an embodiment of the invention. 
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[22] Figure 7A illustrates an example for the reconstruction of a three-dimensional source 
distribution, in which the y=0 slice of the reconstructed helix (for p=l/2 with the 
appropriate filter) is shown in accordance with an embodiment of the invention. 

[23] Figure 7B illustrates an example of the reconstruction of a three-dimensional source 
distribution, in which the y=0 slice of the reconstructed helix (for p=l/2 with the no filter) 
is shown in accordance with an embodiment of the invention. 

[24] Figure 8 illustrates projections of the reconstructed helix, as shown in Figures 7A and 7B, 
as the imaging volume is rotated by 90 degrees around the x axis in accordance with an 
embodiment of the invention. 

[25] Figure 9 shows an apparatus for reconstructing a three-dimensional image representing a 
source distribution in accordance with an embodiment of the invention. 

[26] Figure 10 shows a detector (Dl) flux in accordance with an embodiment of the invention. 

[27] Figure 11 shows a detector (Dl) interaction probability in accordance with an 
embodiment of the invention. 

DETAILED DESCRIPTION OF THE INVENTION 

[28] The following terminology and notation is referenced in the following discussion. 

> BP backprojection operator 

> FBP filtered backprojection 

> FT Fourier transformation 

> LHS left hand side of an equation 

> RHS right hand side of an equation 

> voxel The smallest distinguishable box-shaped part of a three-dimensional space. 
A particular voxel will be identified by the x, y and z coordinates of one of its eight 
corners, or perhaps its centre. The term is generally used in three dimensional 
modeling and sometimes generalized to higher dimensions. 
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> SO Dirac delta function. The function is defined by the relation 

f (0) = f dxf(x) 6(x) 

J for all integrable functions f(x). An alternative 

6 fx) = lim -pi-exp[ - 4 I 
definitionis U ^° «^ *\ ^ I . 

> e x , e y , e z unit vectors corresponding to the x, y, and z coordinates, respectively. 



[29] A mathematical model is proposed for Compton cameras. This model describes the 
expected data from a Compton camera in terms of the source distribution of 
radiopharmaceutical in the patient. In this model, a parameter p is introduced that 
describes the effect of the distance between the source and the primary (scatter) detector. 
A detailed justification of this model is given in Appendix A. The analysis in Appendix 
A indicates that the p=0 case corresponds most closely to actual Compton cameras. 
Backprojection operators (with a distance-dependent index q that is analogous to the 
index p of the model) are then introduced. The analysis of these backprojection operators 
requires the evaluation of an integral that is discussed in Appendix B. Appendix C 
demonstrates how one can select the index q to produce delta functions for p=l/2 and 3/2, 
thereby simplifying the task of the inversion and eliminating the complicated summations 
over spherical harmonics. The backprojection operators are combined with filters (that 
are dependent on camera geometry) and analytic inversions of the Compton camera 
models are found. Because the case p=0 is not solved, one must consider these analytic 
inversions as approximations. The numerical implementation and noise immunity of the 
algorithms are discussed and exemplary results are presented. The implications of these 
algorithms for the design of Compton camera geometry are discussed. 

[30] Figure 1 shows a basic Compton cones geometry in accordance with an embodiment of 
the invention. The basic idea of electronic collimation is illustrated in Figure 1. A 
gamma ray of energy E is emitted at source position (s) 101 (within the patient) and 
propagates to the primary detector (Dl) 105 where it undergoes a Compton scattering that 
deposits energy Ej at position (PI) 111. The scattered gamma-ray propagates to point 
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(P2) 107 in the secondary detector (D2) 109 where it is absorbed by a photoelectric 
interaction and deposits the remaining energy E 2 . (Generally, a detector with low atomic 
number Z is preferable for primary detector 105 so that Compton scatter is more likely; 
whereas, a detector with high Z is preferable for secondary detector 109 so that 
photoelectric absorption is more probable. In addition, primary detector 105 must be 
sufficiently thick for reasonable probability of Compton interaction, but not so thick that 
the scattered gamma ray is reabsorbed in primary detector 105.) In such a Compton 
camera 100, an event is registered whenever coincidence counts are observed in the 
primary and secondary detectors 105,109 that satisfy 

E = E, + E 2 .(1) 

[In reality, the energies can be measured with only limited accuracy, so that events are 
accepted if the sum ( E 1 + E 2 ) is within a window surrounding E.] Thus, for each event, 
Compton camera 100 provides four parameters: the positions of interaction (PI) 111 and 
(P2) 107, and the energies E { and E 2 . Simple relativistic kinematics yields the Compton 
relation 

where m « is the electron mass (51 1 keV) and 8 is the angle of the Compton scattering. 
Two important conclusions can be drawn from equation (2). First, because 
-1 2; cos 8^1, one finds that 

Ei * t — r and , meE r <: E 2 (3) 

1 (2E+m e ) (2E + m e ) 2 V ' 

or, equivalently, \J&\ + 2m e E 1 - E x £ 2E 2 . (4) 

Inequality (4) provides an additional constraint on the acceptance of events, even in cases 
for which the incident energy E is unknown. Second, one observes that equation (2) 
provides an explicit relationship between the angle of scatter and the observed energies 
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a s cos 8 = 



(E 1 + E 2 )E 2 



(5) 



[33] 



[34] 



Compton camera 100 is based on this relationship; namely, the energies Ej and E 2 
determine the angle of scattering (8) 103. Thus, knowledge of the energies Ej and E 2 
implies that the incident gamma ray must have originated from one of the directions on 
the cone defined by a= cos 8; so that, the source (s) 101 must lie on a cone 113 (as 

shown in Figure 1). The basic mathematical problem posed by Compton camera 100 is 
that individual events are associated with cone 113 and not a unique direction (as in 
SPECT or PET). Reconstruction of the source distribution from such information has 
proven significantly more difficult than the standard FBP algorithms for SPECT and PET. 

In addition to relativistic kinematics, further details of the Compton scattering are 
required for the reconstruction of the source distribution from the data of Compton 
camera 100. The probability of scattering at angle (0) 103 is proportional to the 
differential scattering cross section. For a free electron, the Compton scattering cross- 
section is given by the Klein-Nishina formula: 



where r e is the classical electron radius. In the following analysis, the weighting function 
W defined by 



d<*KN(9,E) 
dQ 




(6) 




(7) 



[i + JL^.cose)] 3 
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will be used. The function W gives the relative probability of a photon of energy E 
scattering at angle (6) 103 and has the units of steradian Unfortunately, the electrons in 
detector (Dl) 105 (from which the gamma-ray scatters) are bound, not free. The bound 
state of the atomic electrons causes two effects that influence the performance of 
Compton camera 100. First, the electrons oscillate within the atoms of the detector 
leading to Doppler shifts that affect the relativistic kinematics of equation (5). This 
Doppler shift is characterized by momentum distributions [J(pz)] associated with the 
various atoms in the detector. This effect has been studied and the resultant errors in the 
scatter angle characterized. For low-energy gamma-rays, these errors can be significant. 
For this reason, materials that minimize the Doppler shift should be preferred for detector 
(Dl) 105. Second, scattering at small angles (8) 103. is suppressed if the energy 
transferred to the recoiling electron is insufficient for ejection from the atom. The scatter 
event is forbidden if the excited electron is driven into a previously occupied state, so that 
the cross section becomes dependent on the occupation of outer electron shells of the 
atom (or the Fermi surface of a metal). Extensive tables of the incoherent scattering 
function (which characterizes this suppression for the small angle scattering) are available 
and can be incorporated into the weighting function W. In general, both these effects are 
minor if the initial energy E is large (E>m e ). Despite potential problems for low energies 
(e.g., 140 keV emissions of Tc-99m), both effects are ignored in an embodiment of the 
invention. 

The Compton camera geometry shown in Figure 1 with a small primary detector (Dl) 105 
located in front of a large secondary detector (D2) 109 is a realistic model of the 
prototype geometry built by most researchers. The major limitation is the size of the Dl 
detector 105 which is generally a small, pixelated solid-state detector. A standard 
gamma-ray camera is often used for the D2 detector 109. Since the first development of 
Compton cameras, researchers have realized that this geometry is insufficient for 3D 
reconstruction. The problem can be solved by scanning the camera around the patient, 
thereby effectively increasing the area of the Dl detector 105. Experimentalists tend to 
rotate the camera around the patient in a circular orbit in emulation of SPECT. Some 
theoreticians extend the Dl and D2 detectors 105,109 to infinity and discuss planar 
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detectors. Others surround the patient with spherical or cylindrical banks of detectors. 
Still others rotate infinite planar detectors around the patient. Because each geometry 
potentially requires a different reconstruction algorithm, one might expect a plethora of 
algorithms. However, the filtered backprojection algorithms, in an embodiment of the 
invention, display a significant invariance. The backprojection operation is unchanged by 
the detector geometry; whereas, the filter used on the backprojections depends on the 
geometry of primary detector (Dl) 105 and its scanning trajectory. Consequently, one 
can make a list of relevant geometries (or scanning motions) for primary detector (Dl) 
105 and associate a filter with each geometry. This strategy provides a unified structure 
for development and evaluation of algorithms based on the geometry of primary detector 
(Dl) 105. Five of these geometries are illustrated in Figures 2A-2E. 

[36] Figure 2A shows a rectangular detector in accordance with an embodiment of the 
invention. Figure 2B shows a square prism detector in accordance with an embodiment of 
the invention. Figure 2C shows a cylindrical geometry in accordance with an embodiment 
of the invention. Figure 2D shows a circular geometry in accordance with an embodiment 
of the invention. Figure 2E shows a spherical detector in accordance with an embodiment 
of the invention. Although all of the geometries in Figures 2A-2E have been discussed in 
the literature, few are currently considered practical. Figure 2A shows one of the more 
practical designs, namely, a rectangular Dl detector 201 with dimensions LxA. Figure 2B 
shows a (4-sided) rectilinear box detector 203 with two open ends and each side Lx(2Rc). 
Figure 2C is one of the favorites in the literature: Dl is a cylinder with radius Rc and 
axial length L. Figure 2D shows circular Dl detector 207 with radius L. The limit L-*» 

of Figure 2D corresponds to the infinite planar detector that is popular with theorists. 
Finally, Figure 2E shows a spherical detector consisting of concentric spheres for 
detectors Dl and D2 209,21 1. 

[37] Possible Dl geometries are restricted by the size and shape of the patient. For this work, 
the patient (and, therefore, the source of radiation) will be restricted to a truncated 
cylinder of radius B and axial length A. One observes that B<A for human subjects. For 
some applications (i.e. cardiac or brain imaging), the axial length A may be truncated (to 
the chest or head region). The Dl detector is characterized by a radius Rc- The Dl 
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detector is restricted to regions outside a cylinder of radius Rc>B for all geometries 
except spherical. For spherical detectors, the Dl radius, Rc, must be larger than 

*J A 2 + B 2 where A is the full height of the patient (so that the full patient can fit into the 
spherical detector). 

[38] The secondary D2 detectors are also subject to an important requirement. Not all scatter 
events in the Dl detector will enter the D2 detector. In particular, gamma-rays that are 
scattered back toward the source will miss the D2 detector. Obviously, one wants a large 
D2 detector (as shown in Figure 1) so that as many events as possible will be observed. 
As will be discussed, the equations describing Compton camera performance imply an 
important relation between the Dl and D2 detectors. At each point z£Dl, the D2 

detector is visible (unobstructed by any intervening material) in a solid angle Q ( z , D2 ). 

If this solid angle fi(2,D2) is a hemisphere for every point zGDl, then the 

information is sufficient for the filtered backprojection algorithms in an embodiment of 
the invention. The specific geometry of the D2 detector is not important. The only 
restriction on the D2 geometry is that, from every point in the Dl detector, the D2 
detector must be capable of recording events from a complete hemisphere of scatter 
directions. In Figure 1 , this requirement is satisfied only if the planar D2 detector extends 
to infinity. The algorithms in an embodiment of the invention can still be implemented if 
this condition on D2 geometry is violated, but the resulting reconstructions may be 
suboptimal. 

[39] Figure 3 illustrates three equivalent D2 geometries 301,303,305 that satisfy the 
requirement for a small rectangular Dl detector 351. All the detector geometries in 
Figure 2 A-2E satisfy this requirement. 

[40] Figures 4A and 4B illustrate that the Dl detector geometries in Figures 2 A and 2C can be 
accomplished without the construction of special detector arrays. The same geometry can 
be achieved by scanning with a smaller detector. Furthermore, the geometry of rotating a 
small Ax A detector around the patient is equivalent to the cylindrical geometry of Figure 
2C with L=A. 
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[41] An embodiment of the invention provides a set of FBP algorithms that predict the 
distribution of radiopharmaceutical source distribution within a patient from the data 
acquired by Compton camera 100. For mathematical simplicity, the emissions of the 
radiopharmaceutical source distribution will be assumed monoenergetic with energy E. 
[This assumption can be dropped without affecting the final imaging equation (11). 
However, the notation becomes more complicated.] The distribution of 
radiopharmaceuticals is described by a function f ( x ) that denotes the density of 
radiopharmaceuticals in the patient. The function f has dimensions cm -3 

(emissions/cm^) and is a bounded, positive function with support within a cylinder of 
radius B and axial length A (i.e., the radiopharmaceuticals are contained within the body 
of the patient). 

[42] The data generated by Compton camera 100 are a collection of N coincidence events 
from the primary and secondary detectors 105,109. For SPECT studies, one typically 
collects N « 10 7 counts. Whereas, for Compton camera 100, one expects (at least) a one 
hundred-fold increase in sensitivity, so that N«10 9 . Each event provides four 

parameters that describe the coincident detection: the position z of the Compton 
scattering in primary (Dl) detector 105, the position w of the photoelectric absorption in 
secondary (D2) detector 109, and the energies E { and E 2 deposited in these detectors. In 

summary, a coincident event provides two positions and energies: (z,E l ;w,E 2 ) where 

z E Dl and wED2. For purposes of data analysis, these data will be transformed into a 

new set of parameters: 



(z,E 1 ;w,E 2 )=>(z,n,a f E) 



by the algebraic relations 




, E = E x + E 2 , and 
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z represents the position of Compton scattering in the primary detector, n denotes the 
direction of the y-ray after scattering, a corresponds to the cosine of the angle between 
incident and scattered directions of the y-ray, and E corresponds to the energy of the 
incident y-ray. Each event is characterized by four parameters ( z , n , o , E ), such that 

zGDl , |n|=l , and -1<o<1. (9) 

[43] Thus, the two forms of the data have the same number of dimensions (6) and the 
parameters E x and E 2 are interchangeable with a and E. However, the vectors w and 3 
are subtly related by the geometry of the detector system; \56D2 whereas | n | = 1 . The 
vector n denotes the direction of the y-ray propagation after scatter. At each point 
zEDl some directions n will not intersect the detector D2 109 and, therefore, are not 

sampled by the detector system. For this reason, Compton camera 100 is often called a 
"limited angle" imaging system and the data are called "incomplete." From this 
argument, many researchers have concluded that backprojection algorithms are not 
suitable for Compton cameras and that alternative (e.g., maximum likelihood) algorithms 
are the only appropriate method for reconstruction. One of the important ideas, in an 
embodiment of the invention, is that the data are complete if the vector ii is sampled on a 
hemisphere at each point zEDl [see the discussion following equation (12)]. Many 
current Compton camera designs meet this criterion. 

[44] [If the energy E of the incident y-ray is known, the energy detected in secondary detector 
D2 109 becomes redundant. As a result, one could accept scatter events in detector D2 
109, rather than requiring total absorption in a photoelectric interaction. This strategy is 
often proposed for high-energy radiation because the probability of photoelectric 
absorption in the secondary detector drops for higher energies. However, higher rates of 
accidental coincidence counts may result. In an embodiment of the invention, the 
reconstruction algorithms require the cosine of the scattering angle and the energy of the 
incident y-ray E. The method by which these parameters are determined does not affect 
the algorithms.] 
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[45] The list of events from the camera is next combined into a histogram of events 
g ( z , n , a , E ) that describes the imaging process. Because n and a are expressed in 
dimensionless parameters (steradians and cosines, respectively), the resulting histogram g 

-2 

function has dimensions cm corresponding to the density of observed events at the 
position z on the 2D surface Dl . The observed function g(2,n,a,E)is related to the 
source distribution f ( £ ) by the integral equation: 



g(z,n,a,E) = W(a.E) f f j 

= W(o,E) ///d 3 X f(X) 



2jtRP| 2-x | 2 " p 



2jtRP| z-x | ! - p 



6 



n • ( z - x ) 
I | 



-a 



6 [ n • ( z - x )- a| z-x | ] 



(10) 



[46] The physical justification and limitations of equation (10) are discussed in Appendix A. 
For all realistic imaging processes, the parameter p should vanish; the case p=0 implies 

-2 

that the radiation dies off as r as Dl detector 105 moves away from the source. For 
technical reasons, however, analytic inversion for the p=0 case (and more generally any 
integer value of p) is difficult. On the other hand, the p = ±1/2 and 3/2 cases yield 
analytic solutions. The imaging equation (10) is a generalization that includes these 
analytic inversions. The function W is the dimensionless weighting function [defined in 
equation (7)] that accounts for the different count rates at various scattering angles. One 
can immediately absorb the function W into the LHS of equation (10) and write that 



W(o.E) l)) dX HXJ 2 3C RP|2-X| 2 -" 



6 



n-(2-S) 
12-21 



- o 



= fffdH f(X) 



2nRP| i-a | l_p 



6[n-(z-x)-o| z-x |] (11) 
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[47] All energy dependence in equation (10) arises from the function W; so that, the resulting 
ratio g/W is independent of the total energy E. Equation (11) will be called the Compton 
imaging equation with index p and C p is the Compton camera operator. This notation is 

slightly misleading because, once the function W is factored out, equation (11) actually 
describes an isotropic scattering process. However, the important idea, namely that the 
angle of scatter is determined by the energies, is retained. Straightforward dimensional 

-3 -2 

analysis confirms that f has units of cm and g has units cm . One should also note the 
symmetry 

g(z,rl,a,E) _ g(z,-n,-a,E) 
W(o,E) W(-o,E) 

[481 This symmetry is very important because, at any point z, in the primary detector, half of 
the directions n are generally excluded from measurement. Equation (12) guarantees that 
the remaining directions n are sufficient for the determination of g over the full domain 
of n. Thus, equation (12) is the basis of the geometric condition that detector (D2) 109 
must cover a hemisphere of directions n for each point in detector (Dl) 105. Another 
important property of equation (11) is that the integral on the RHS is a convolution 
integral in ( z - x ). Normally, convolution integrals are easily inverted with a Fourier 
transformation (FT). Unfortunately, the function g is known only on detector (Dl) 105, 
not over the entire 3D range of z and, therefore, cannot be Fourier transformed. Thus, 
the most obvious inversion technique is not feasible. 

[49] Figure 5 shows a conceptual flowchart for filtered backprojection (FBP) algorithms in 
accordance with an embodiment of the invention. The inverse problem for a Compton 
camera [of index p] can now be simply stated: one must find the function f in equation 
(11) that generates the observed function g. The inversion strategy employed, in an 
embodiment of the invention, is the standard backprojection strategy as illustrated in 
Figure 5. The source distribution f(X) 501 produces the Compton camera data 

g(z,n,a,E) 503 by experimental observation [as modeled by equation (10) with 
index p]. The inversion requires the selection of an appropriate combination of weighted 
backprojection operator (BP) 505 and filter function (O) 507. An important step in this 
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process is the selection of the correct backprojection operator. For the success of this 
inversion strategy, the FT of the backprojected data must be proportional to the FT of the 
original source distribution. The relation between these FTs is an integral equation 
[which is given in equation (24)] and the strategy can only succeed if the kernel in this 
integral equation is proportional to a delta function. The introduction of the apparently ad 
hoc weighting indices p [in equation (10)] and q [in equation (13) below] is required 
because the correct combination of p and q yields a delta function in equation (24) that 
renders filtered backprojection feasible. 

In an embodiment of the invention, the backprojection operator (with index q) of g is 
defined by the equation: 

^q[w] = b q( § ) 

(13) 



f d 2 2 ( d2fi \ da g( 2 ; R '°f) _i b \*dlzll 

Jdi J { x W(a,E) 23tR^|z-S| 2 - <1 



|z-s| -° 



where Dl corresponds to the surface of primary detector 105. The integrals in equation 
(13) are analogous to the standard backprojection done in SPECT and PET; except that, 
rather than backprojecting the data onto a line, one backprojects onto a cone (described 
by the delta function). [Mathematically, the operator BP p is the adjoint of C p .] The 
index q is analogous to the index p in the imaging equation (11) and permits weighting of 
the backprojection with various powers of | s - z | as one projects the data inward from 

the point of scattering z in detector Dl 105. The motivation for the introduction of the 
index q will become apparent later in the calculation [see equations (29) and (30) and 
Appendix C]. The domain ( s ) of the functions b q is the 3D region corresponding to the 

possible support of f and the functions b q have units of cm 3 for all values of q. Thus, all 

the functions b q are candidates for the inversion of Compton camera imaging operator 

C p . The integral over n in equation (13) requires integration over the entire 2 -sphere. As 

previously noted, the function g is seldom observed for ii over the full 2-sphere. 
However, if the function g is observed over a hemisphere, one can use the symmetry of 
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equation (12) for evaluation of g in the unobservable directions. In practical terms, this 
means that the integral over n in equation (13) can be evaluated using an integral over the 
"observed" directions n on a hemisphere; the contribution to equation (13) from the 
"unobserved" directions merely doubles the value of the resultant integral. 

[51] The backprojection operators will be examined using the functions h q „ defined by 



h qp (s) = BP q C p [f] 

.///A f(x)/ Di d- Z 



(2*) 2 R»^ + P|2-§| 2 -< , |z-2| 2 - p 



(14) 



fB|-l -1 



da 6 



n • ( z - s ) 
12-31 



- a 



r -» 



n * ( z - x ) 
z-X I 



- a 



[52] One might naively hope that h q p = f for some p and q because, in that case, BP q would 

be the inverse operator for C p . Despite the fact that such hopes are not met, the 

evaluation of the functions h q p will lead to useful FBP algorithms in the same manner 

that filtered backprojection arises in SPECT and PET. Performing the integral over o in 
equation (14), one finds 



h,.p(§)=///^ f(x) / d /z 



(23i) 2 R^ + P| 3-2 | 2 ^| x-z | 2 " p 



(15) 



f d 2 n 
|B|=1 



n • ( z- s ) fi • ( z - x ) 
I s - z I I x- 2 I 



[53] Next, the integral over the 2-sphere n can be evaluated by a simple calculation 
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Zx n 



2n 1 



f d 2 n 6(n-v) = f d<p jd0 sin 9 6(|v|cos6) = J dtp f d\i 6(|v|n) 
|a|=i oo o-i 

= 2 3 t|dn6(|^| l x) = ^|dn6( l i) = ^ 

so that in equation (15) one finds 

h q.p( § )=/// d3 * f (^) 



(16) 



L 



(17) 



d 2 z 



27tR^ + P| s-z f^l x-z | 2 " p 



(£-2) j3-2)_ 
lx-21 |s-2| 



[54] Next, the 3D Fourier transformations (F and Hq s p) of the functions f and h q p are defined 
by the equations 



and 



F(v) = /J/d 3 x f(x) exp[-2ra(v-x)] 
H q.p( 3 ) = /// d3§ h qp( § ) exp[-23ci(63-s)] 



(18) 
(19) 



[55] Thus, equation (1 7) can be written in the spatial-frequency domain as 

H q p (63) = jfj d 3 tf F(v) jff d 3 x JJJd 3 s exp [ 2ra v • X - 2jii 3 • s ] 

I 



(20) 



d 2 z 



2jiR^ + p| s-z | 2 ^| x-z | 2 " p 

[56] One can translate the dummy variables of integration s and x, 
X -* x + 2 and s -» s + t , 
so that 



(*-*) 

|x-z| | s-z I 
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H qp (<3) = fffd 3 v¥{1) J^d 2 ? exp[2ati(v-S)-z] 

///* /// 



(21) 



. exp [ 2ni v • x - 2ni 63 • s 1 
d 3 s ^ 1 



% L 

x I 15 



[57] Thus, one finds that the integrand on the RHS of equation (21) contains two independent 
terms: 



H,. P (*.S) =///d 3 * ///d 3 s 



exp [ 2jii v • x - 2ni 65 • s ] 



and 

such that 



2*R^ + P| s | 2 ^| x | 2 ~ p 
D(v) = J Dl d2 2 exp(-2jiiv-z), 



1*1 1*1 



(22) 



(23) 



(24) 



[58] The function M qp contains all the important information about the backprojection 

operation, whereas the function D is a simple Fourier transform over the surface of 
primary detector Dl 105. In the subsequent analysis, one finds that the properties of 
M q p determine the weighting index q of the backprojection operator and the function D 

determines the appropriate filter function. The function M q p will be examined first. 

[59] The function M q p is evaluated by explicit integration in Appendix B. Unfortunately, the 
integral in equation (22) is only convergent if -1 <p< 1 and -1 <q< 1. The lower 
limits on convergence (-l<p and -l<q) arise from the small distance behavior of the 
integrand ( | x | -> 0 and | 3 | 0 ) and are intrinsic to the problem. However, the upper 
limits (p<l and q<l) arise from the large distance behavior of the integrand 
( | x | — > oo and | s | -* oo j. For all practical applications, both the source distribution f 

and the backprojections h q p have compact support, so that the behavior at extremely 
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large distances should be unimportant. This suggests that the integral can be regularized ' 
and the divergences at large distances suppressed. In particular, one can choose 



which is convergent for -1< p and -1< q. For all positive values of Qthe additional term 



Furthermore, as shown in Appendix B, the limit as e -» 0 is well-defined. The 
regularized integral in equation (25) will henceforth be used as the definition of M q p . 

This regularization is significant because the cases q=3/2 and 5/2 will be important in the 
subsequent inversions as will be discussed. 

[60] For mathematical convenience, a unit vector on the 2-sphere will be denoted by 

j~r = C2k = Q ( e k • <Pk ) = cos 0 k e 3 + sin 6 k sin qp k e 2 + sin 6 k cos q) k e { . (26) 
| k | 

[61] The evaluation of the function M q p in Appendix B yields that 



M q>p (v,3) = Hm o J// d 3x /// 




(25) 



in the exponential eliminates the divergences in the integral as | % \ —■ 



co and | s | -* 



M qp (v.S) = 




(27) 





(28) 



one observes (Appendix C) that 
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6 2 (^v-^o) 

2^R||65| 3/2 |>?| 



and M 1/2i3/2 (v,o)) = — — — — ^ . (30) 



[62] The delta functions in equations (29) and (30) are extremely suggestive. The inversion 
algorithms for the Compton camera equation with p=l/2 and 3/2 are a direct result of 
these equations. Indeed, the discovery of functions M qp (or combinations of such 

functions) that produce delta functions of the form 6 2 ( - ) is an important idea in 
an embodiment of the invention. 

[63] For the case p = -1/2, the functions M 5/2 _ 1/2 and M, /2 _ 1/2 prove useful. In particular, 
one finds that 

-siijij»j*F £|£Hfj x,v[a.) [v?(A.)r 

(3D 
and 

(32) 

so that 

2 / — • — » > 

This last identity provides the basis of inversion for the p = -1/2 case. 

[64] The three equations (29), (30) and (33) (with various parameters p and q) all yield the 
same basic expression that will be used in the subsequent inversions. For convenience, 
that mathematical expression will be defined as 
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6 2 (fr,-3») 



[65] For those combinations of p and q that produce jOt k , one can rewrite equation (24) as 

ft k (<5) = |J/d 3 v F(v) M k (v,i5) D(<3-v) (35) 

OO 

sothat ^R^H^AJojQJ = f dv |v| k f(v Q w ) D[(to-v)Q to ] (36) 

OO 

[66] [In equation (36) the identity 



OO 



rrjd 3 v= fdv|v| 2 frd 2 o v = ijdv|vi 2 jjd 2 Q v on 

is used so that the range of the integral is (-oo, oo) rather than [0, oo).] An important 
observation is that the RHS of equation (36) is a ID convolution integral. Thus, standard 
Fourier techniques permit the analytic inversion of equation (36) and the recovery of the 

function F from the backprojected data,H k . The function D is obviously important in 
such an inversion. 

[67] The function D is determined by the geometry of the primary detector Dl 105. The 
definition of D [equation (25)] suggests that, for many geometries, D will contain delta 
functions. Such delta functions may simplify the integral in equation (36). However, not 
all geometries will produce delta functions and care is required in many geometries that 
do produce delta functions. 

[68] In an embodiment of the invention, five basic detector geometries are examined (Figures 
2A-2E). In Table 1, each of these detector geometries is described and the function D 
displayed. 
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Geometry 


Dl Surface 


D(v) 


Jim D(v] 

(L»4max(A,R c )j 


1. Rectangular Detector 


-L/2<z<U2.-A/2sysA/2 ,1 = 1^ 


X «p(2raR t ? ? I ) 


D,(*) = Afij^JsncfsL*-^) 


2. Square Prism Detector 


-Ul<i<\J2 , -R c sysR c , x=±R Q 
-V2<z<V2 , -R c sxsR c , y=*R c 


D 2 (v) = 4R e Liittc(nL?-« 1 ) 

lfitt(2*MI,)<»(2»M«i) 

■ 

.ic(l«,ii 1 |«.(i«,!y 


«*(2iM? T )«(*MM 


3. Cylindrical Detector 


-L/2<z<L/2 .x 2 + y 2 = R 2 


D J («) = 2«R c Laiie(«Lil? 1 ) 




4. Circular Detectors 


Os(y 2 + z 2 )sL 2 .xrRc 




D 4 (*) = 4L6(^*(*tf) 
xeipj^M-*,) 


5. Spherical Detector 


x 2 + y 2 +z 2 =R 2 


D 5 (5) = 4^^(2^(51) 


D 5 (i^) = 4^^(2^191) 



TABLE 1 : The function D(v) for different camera geometries 
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[69] Detector geometry 1 (first row of Table 1) is a rectangular area with dimensions L x □. 
The actual detector need not have physical dimensions L x A, one could scan a much 
smaller detector over the region (Figure 4A) and the resulting function D would be the 
same. Detector geometry 2 (second row of Table 1) is a square prism (an elongated box 
without top or bottom). This geometry is a simple application of the first geometry; the 
function D can be evaluated by a summation of the terms for the four sides. Detector 
geometry 3 (third row of Table 1) is a cylinder of radius Rc and length L. Once again, the 
primary detector Dl 105 need not be a cylinder if the surface Dl is sampled by scanning. 
One might build a ring detector and scan it along the axis or one might build a narrow 
strip detector 403 and rotate it around the axis (as shown in Figure 4B). Either 
configuration would yield the same function D. Detector geometry 4 (fourth row of Table 
1) is a large circular detector of radius L. Finally, detector geometry 5 (bottom row of 
Table 1) is a spherical detector with radius Rc surrounding the imaging region. 

[70] Detector geometries 1 through 4 will be examined in the limit as L-*oo. Recalling that 

Jim o [Lsinc(3tLv)] = 6(v) (38) 
J 2 (2jcLv) 



and lim 

L oo 



* L 2jcLv 



= 6(v) (39) 

one observes that these limits yield the delta functions needed for the simple inversion of 
equation (36). In reality, of course, no detector can be extended to infinity. However, 
because the region of source reconstruction has finite dimensions (cylindrical radius Rc 
and axial length A where generally Rc<A), one need not sample the spatial frequency in 
steps smaller than (Av) « V( 2A ). The width of the sine function in equation (38) is 
approximately (Av) « 2/L. Thus, if L > 4A, then L should be sufficiently large that the 

function D can be approximated as a delta function. [A similar argument applies for the 
Bessel function in equation (39).] In particular, the last column of Table 1 displays the 
L-*oo limit of the function D for detector geometries 1-4. For the spherical detector 
geometry (5), no limit yields a delta function. 
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[71] Equation (36) provides the prototypical equation for all the inversion algorithms in an 
embodiment of the invention. Using the delta function approximations for Di of detector 
geometries 1 through 3, one finds that 

d,[(— vjn.] = k, »[(—,)(*..«,)]- K ;°< m ;;> m 

I • e E I 

where the constants Kj are given by 

K, = A , K 2 = 8R C ,and K 3 = 2jiR c . (41) 
[72] Substitution of equation (40) into equation (36) yields 

i ik 

4*R3|(o| 2+k ft k (<o£2 u> ) = , FfcoflL) (42) 

I Q » * S Z I 

or F (v) = ^|v||v-e z |ft k (v) (43) 

[73] Equation (43) provides a filtered backprojection algorithm for the source distribution F in 
detector geometries 1, 2, and 3. [The second equality in equation (40) presumes that 

I ^id'^z |*°- Consequently, the results in equation (43) are not rigorous on the 
| v ■ e z | = 0 plane. Nonetheless, equation (43) is satisfied for v arbitrarily close to the 
| v • e 2 | = 0 plane. Thus, the continuity of the function F suggests that equation (43) will 
also be valid on the | v • e z | = 0 plane.] Surprisingly, the algorithm does not depend on 
the index k because the same power appears on both sides of equation (42) and is 
canceled. For circular detector geometry 4, the delta function approximation yields a 
similar result, namely, 

F(v) = | v |v^-ay) 2 + (*-8.) 2 Ak(v) (44) 

[74] Equations (43) and (44) demonstrate that, if the detector geometry produces a function D 
containing a delta function, a simple filtered backprojection algorithm follows 
immediately. 
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[75] Unfortunately, the spherical detector geometry (5) does not produce a function D5 
containing a delta function. The substitution of D5. (from Table 1) in equation (36) gives 
the 1 D convolution equation 

00 

2R2|a ) | 2+k ft lc (a>3 a ,) = fdv|v| k f(v^) 2R C sine [ 2*R C ( o - v ) ] (45) 

-OO 

[76] Fourier analysis provides a straightforward method for inversion of such a convolution 
equation. However, a property of the sine function provides an even more powerful 
method. If one defines 

s(co) s 2R C sine ( 2jiR c o> ) (46) 

CD 

then s ( v ) = ( do) s ( v-tn ) s ( a) ) (47) 

-o© 

[77] [Equation (47) follows immediately from the fact the Fourier transform of s is equal to 
unity on the interval ( -R c , R c ) and vanishes elsewhere. As a result, the square of the 

Fourier transform of s is equal to itself. By the convolution theorem, therefore, s must 
satisfy equation (47).] Thus, by taking a convolution integral with s on both sides of 
equation (45), one finds that 

00 00 

I da>2R2|co| 2+k ft k ((o^ J) )s(v-co) = f dco | a> | k F ( a> ^ ) s(v-o) (48) 

-OO -00 

[78] The obvious solution of equation (48) is 

F(S) = 2R2|(o| 2 ft k (aS) (49) 
[79] Of course, solutions of the homogeneous integral equation 

OO 

|dco|a>| k F^oQj s(v-cd) =0 (50) 

-00 

can be added to equation (49). The selection of a particular homogeneous solution 
requires an additional criterion. For all practical purposes, the homogeneous term can be 
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assumed to vanish and equation (49) gives the resulting solution. Equation (48) implies a 
consistency condition on the function H k ; namely, 

00 

|v| 2+k ft k (vQ l0 ) = f dcoH^ftkfoqJ 2R c sinc[23tR c (v-o>)] . (51) 

—00 

[80] Filtered backprojection algorithms can now be stated for the Compton imaging equation 
of indices p = ±1/2 and 3/2. These FBP algorithms are based on five assumptions that 
have been stated earlier, are implicit in the previous analysis, and are repeated here for 
completeness. First, the data from the Compton camera take the form of a function 
g(z,n,o,E) that is known on the domain defined in equation (9). Second, the 
primary detector surface Dl is assumed to be one of the five geometries described in 
Table 1 (Figures 2A-2E). Third, the source distribution f ( x ) has support only within a 

compact region. For detector geometries 1 -4, the support of f is a cylindrical region of 
radius Rc and axial length A. For spherical detector geometry (5), the support of f is a 
spherical region of radius Rc- Fourth, the parameter L in geometries 1 through 4 must 
satisfy the condition L» 4 max ( A , R c ). Finally, the functions g and f are related by 
the Compton imaging equation (11) with index p = ±1/2 or 3/2; that is 

g(z > n,o,E) 

W(o.E) pL J 1 ' 

[81] Of course, the FBP algorithms will be different for the various indices p and detector 
geometries. All the FBP algorithms can be divided into four steps. The first step 
involves the backprojection of the data and is dependent only on the index p. The second 
step requires the 3D Fourier transformation of the backprojected images. The third step 
involves the filtering of the backprojection data and depends only on the detector 
geometry. The final step is a simple 3D Fourier transformation that recovers the source 
distribution f and is independent of either index p or detector geometry. 

[82] Figure 6 shows a detailed flowchart of filtered backprojection (FBP) algorithms in 
accordance with an embodiment of the invention. Figure 6 delineates the steps (steps 
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603, 605, 607, and 609) that constitute FBP algorithms in accordance with an 
embodiment of the invention. 

STEP 1 (step 603): Backprojections for index p 

[83] In step 603, the Compton camera data g ( 2 , n , o , E ) (representing input data 601 and 
generated by Compton imaging equation of index p) are mapped by backprojection into 
various functions b q (x) = BP q [^w] [equation (13)]. For each index p, different 
indices q are required for the backprojected functions b q (x). In particular, the 
algorithms for the three cases p = ±1/2 and 3/2 require the following backprojections: 

P=3/2 requires b l/2 (x) [ = h 1/2f 3/2 ( x ) ] , 

E^J/2 requires b 3/2 (x) [ = h 3/2i 1/2 (x)], 
and 

P = -1/2 requires bothb 5/2 ( x ) [ = h 5/2| _ 1/2 ( x ) ] 

and b 1/2 (2) [ = hi/2.-i«(*)]- 

STEP 2 (step 605): 3D Fourier Transformation of the Backprojected Data 

[84] In step 605, the backprojected data are Fourier transformed and produce an intermediate 
function H(v). For notational convenience, the 3D Fourier transform of the 
backprojected data is defined by 

B q (v) = ///d 3 X b q (x) exp(-2jiiv-x) (52) 
[85] Algorithms for the three cases p = ±1/2 and 3/2 are then given by: 

Er3/2 H(v)= B 1/2 (*) [ = H 1/2>3/2 ( v ) ] , (53) 
H(v) = B 3/2 (v) [= H 3/21/2 (v)] , (54) 
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and p = -1/2 



H(v)= Bs, 2 (v) + 



(2*R c |v|) : 




(55) 



= H 5/2 _ l/2 (v) + 



(2^R c |v|) : 



2 Hia.-iaO • 



[86] For each value of p, the function H corresponds to the combination of functions Hq 5 p that 



produced M k in equations (29), (30) and (33). The function H, therefore, satisfies the 
integral equation (36) as previously discussed. 

STEP 3 (step 607): Geometric Filtering 

[87] In step 607, the intermediate function H is converted to the Fourier transform of the 
source distribution, F, by the appropriate filtering (as previously discussed). The filtering 
is dependent on the Dl detector geometry. Thus, the function F is given by 



F(v) = <D(v) H(v) 
where the function O ( v ) is a filter function determined by the Dl geometry. 



(56) 



For detector geometries (i=l-3) (corresponding to Figures 2A, 2B, and 2C): 




(57) 



where Ki are given by equations (41). 



For detector geometry (4) (corresponding to Figure 2D): 





(58) 



For the spherical detector geometry (5) (corresponding to Figure 2E\ 
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the function <I> is given by 

= 2R2|v| 2 . (59) 
STEP 4 (step 609): Inverse 3D Fourier Transformation 

[88] In step 609, the source distribution f is recovered from F by inverse Fourier 
transformation [of equation (18)]: 

f (X) = jjj d 3 v F(v) exp[2*i(v-x)] (60) 

[89] Steps 601-609 constitute the FBP algorithms in accordance with an embodiment of the 
invention. In condensed operator notation, these algorithms can be written as 

Weighted 
Filter Backprojection 

= ^ [*(*)] F 3 bp^ 

Cm = Tj l [*{*)] Ts BP3/2 (61) 



<D(v) 



T-x BP, 



1/2 



where T 3 denotes the 3D Fourier transformation and O ( v ) is the filter function 
associated with the Dl detector geometry. 

[90] Throughout the derivation of these algorithms, Compton camera data 601 are tabulated as 
a histogram function g(z,n,o,E). In reality, of course, the data are typically 
accumulated as a list of N individual (coincident) events of the form 

( 2j ,nj , (Jj , Ej) for i = 1 N 

[911 Organizing these events into a histogram g ( z , 5 , a , E ) is one method of handling the 
data. An alternative strategy is typically more practical. If one estimates the function 
g ( z , n . a , E ) by 
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g(z,n,o,E) = i | i 6(z-z i ) ftfa-fij) 6(0-0*) 6(E-E S ), (62) 
the (linear) backprojection operation in Step 1 can be performed on each individual event 
and accumulated into the functions b q ( x ) as 



b q (x) = | i fiP < 



6(z-Zi) 6(n-ni) 6(o-a;) d(E-Ej) 



(63) 



[92] Thus, one jumps directly from the list mode data to the backprojection function b q ( X ). 
The backprojection operations BP q require much more computation than binning the data 

into voxels of a histogram. Consequently, if many events fall into the same voxel of g, 
the accumulation of events into a histogram would reduce the computational burden. 
Assuming that events are limited to a narrow energy window, the domain of the function 
g(z,n,a,E 0 ) is five dimensional, whereas, the domain of b q (x) is three 

dimensional. This difference is important. If one samples the functions g ( 1 , n , o , E 0 ) 
and b q (x) with 128 (2?) subdivisions for each dimension, then g requires 3.4x1 0*0 
(2^5) voxels, whereas, b q requires only 2.1x10^ (2^1) voxels. For a typical acquisition 
with Compton camera 100, the current estimate gives the number of accumulated counts 
as N« 10 7 - 10 9 . If 10 9 counts are distributed among 3.4xl0 10 voxels, it is unlikely 
many events will fall into the same voxel. Thus, any computational advantage of 
rebinning the data into a histogram g is lost. Thus, one concludes that the histogram 
function g(z,n,o,E 0 ) should not be explicitly calculated. Instead, the list mode data 

should be accumulated directly into the appropriate backprojection functions b q ( x ) as 
prescribed by equation (63). [For p=l/2 or 3/2 only one function b q is required; whereas, 
for p = -1/2 two functions b q are required.] This strategy requires significantly less 
computer memory and should not slow the reconstruction. 

[93] A major consideration in image reconstruction is noise immunity. FBP algorithms for 
CT, SPECT, and PET are surprisingly immune to noisy data, but the FBP algorithms for 
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Compton camera 100 have worrisome features that suggest noise immunity could be a 
problem. In the (2D) FBP algorithms for CT\ SPECT, and PET, a ramp filter | v | 
appears in the reconstruction. This ramp filter enhances high spatial-frequencies and 
noise. This ramp filter is generally multiplied by a (low-pass) noise filter {e.g., a 
Butterworth filter) that suppresses the high-frequency noise. In general, the (3D) FBP 

algorithms for Compton cameras contain a quadratic filter | v | 2 rather than a linear 

(ramp) filter and, therefore, enhance the high-frequency noise even more. One might 
expect, therefore, that a stronger noise filter would be required for Compton camera 
reconstruction. Inherently, a noise filter causes a loss of image resolution, so that one 
expects either noisier images or a loss of resolution in the reconstruction. Fortunately, 
there is an alternative method of handling the quadratic filter that apparently avoids this 
problem. Noting that 

§ 3 - , [-4 3t 2|v| 2 ]^=V 2 > (64) 
one can write the filtered backprojection reconstruction algorithms [equation (61)] of the 
previous section as 

O(v) 



£3/2 - 



^1/2 - T: 



<F 3 (-V 2 )BI\ 




(65) 



|(-V 2 )BP^ + ^gP l/2 |. 



The filter function O ( v ]f\ v | 2 is well-behaved at large spatial-frequencies and does not 

amplify noise. A numerical implementation of the Laplacians in equation (65) has been 
used successfully for data containing fewer than a million counts (so that significant 
statistical noise would be expected) without any noise filter. 

Another concern about noise immunity arises from the weighting applied in the 
backprojection process. The weighted backprojection operators BP q add different 
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amounts to the images depending on the distance from point of scattering in Dl detector 
105. If q<2, the backprojection operator amplifies the signal near the point of detection in 
detector (Dl) 105 and suppresses the signal far away. On the other hand, if q>2, the 
backprojection operator suppresses the signal near the point of detection and amplifies it 
far away. The latter case is bothersome because experience with attenuation correction 
algorithms implies that amplification of signals far from the point of detection amplifies 
noise. Only a small region is close to the point of scatter (in Dl), whereas most of the 
universe is far from this point. Thus, amplification of the signal far from the point of 
detection implies amplification almost everywhere. The algorithm must subtly cancel 
such large signals in the far field for successful reconstruction. Noise disrupts this subtle 
cancellation and produces significant artifacts in the reconstruction. Thus, the 

backprojection operator BP 5/2 that appears in the inversion for p = -1/2 ( CZln ) causes 
concern. This potential problem is best evaluated with computer studies that simulate 
statistical noise. 

[96] Because of these concerns, preliminary computer simulations were run to test whether the 
algorithms amplified noise. These initial simulations were performed for the spherical 
detector geometry (5) with the Dl detector radius of 6.3 cm. The reconstruction volume 
was (128x128x128) voxels with 1 mm on each voxel side. Inside this volume, we placed 
a helix source phantom 6 cm long (z axis) and having diameter 1 cm. This phanton can 
be visualized as a simple helical wire defined by the affine parameter t: 

z = 6 t-3 

x = 0.5 cos ( 4;rt ) 

y = 0.5 sin ( 4xt ) (66) 
where 0<t<l. 

[97] The study involved the simulation of 1 million detected scatter events. Each event 
required the tracking of a photon emitted in a random direction from a point on the helix, 
then propagated to the Dl detector 105, scattered at an angle sampled from the function 
W, and then propagated to D2 detector 109 where the detector was assumed thick enough 
to assure absorption. Thus, the simulation produced data appropriate for the physical 
(p=0) model of the Compton camera. The scattering was assumed isotropic (W=l, rather 
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than Klein-Nishina) and no Doppler broadening effects were included. (Later simulations 
using the Klein-Nishina cross-section provided virtually identical results.) The spatial 
and energy resolution of the detectors were assumed perfect, so that the only sources of 
error were the statistical noise and the index p of the reconstruction algorithm. 
Reconstructions were performed for the filtered backprojection algorithms described by 
p=l/2, 3/2, and -1/2. No noise filters were used in the reconstructions. The algorithms 
all produced comparable results -- with the conclusion that p=l/2 gave the marginally 
best results, 3/2 next best, and -1/2 the least successful results. The observed resolutions 
for all these FBP reconstructions were <1.5 mm — essentially limited by the voxel size 
(1mm) of the simulation. Figures 7 and 8 illustrate the results of the FBP reconstruction 
for p=l/2. 

[98] Figures 7 A and 7B show the reconstruction results on a slice through the y=0 plane. 
Figure 7 A shows the results of FBP with the appropriate (Laplacian) filter; the four points 
where the helix intersects the y=0 plane are clearly visible. Figure 7B shows the result of 
backprojection without the filter; the blurred image is comparable to the results of Rohe 
and Valentine. Visualization of the 3D distribution reconstructed from the helix source is 
difficult. 

[99] Figure 8 shows a sequence of projections of the reconstructed distribution from various 
angles as the imaging volume is rotated 90 degrees around the x axis. The helix structure 
is clearly visible. The "bead" structure observed on the helix is a partial volume effect 
caused by superimposing the finite voxel size (1mm) on the infinitely thin helical wire. 
These images suggest that neither the statistical noise nor the artificial index p pose 
significant problems for the FBP algorithms proposed in an embodiment of the invention. 
However, the unphysical assumptions (perfect spatial and energy resolution) leave major 
questions about the application of the algorithms to realistic data. More extensive 
simulations are required and currently underway. 

[100] In an embodiment of the invention, an integral equation (11) was introduced that 
describes the imaging process in Compton cameras. That integral equation was 
generalized with an index p that characterizes the response of the camera to sources at 
different distances from the primary detector. The case p=0 (corresponding to the 
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radiation flux decreasing with distance as r ) is an accurate model of Compton camera 
performance. In an embodiment of the invention, analytic inversions were found for the 
cases p = ±1/2 and 3/2 for certain detector geometries. These analytic inversions take the 

form of FBP algorithms. If the inverse Cq 1 were known, the analytic problem would be 
solved. Because Cq 1 is not known, one must hope that cyj^, Cr} /2 , or C3/2 will serve as 
an acceptable approximation of Cq 1 . (Attempts to find a FBP algorithm for Cq 1 are still 
under study. However, the most promising result yields an infinite series of different 
weighted backprojections and filters. Aside from serious questions of convergence, such 

expressions for £q 1 completely unsuitable for practical applications.) If one uses C' 1 

rather than Cq 1 for image reconstruction, the resulting algorithm will produce systematic 
errors. These errors can be studied either analytically or experimentally with computer 
simulations. Analytically, one examines the difference 

E p = C^Co- I (67) 
where E p f gives the error in the reconstruction of a source distribution f. The analytic 

evaluation of E p is not very difficult — the calculation is basically the same as previously 

discussed. The calculations of Hp,q, Mq, p and D recur and can be evaluated in the same 

manner. A (chi-squared) norm of the error operator (i.e., |E p |^ can be found, but such 

norms are notoriously bad indicators of performance in image reconstruction. 
Experimental computer simulations and observer performance studies are required for 
evaluation of these algorithms. Such studies are now underway. 

[101] In addition to the actual FBP algorithms, the mathematical analysis provides significant 
insights into the geometrical design of Compton cameras. Typically, the primary detectors 
of corresponding Compton cameras are small solid-state detectors. These small detectors 
may be sufficient for proving the concept because the small detector can be moved 
around the patient (or experimental phantom) in simulation of an extended primary (Dl) 
detector. However, most researchers then discuss moving the camera around the patient 
in a circular orbit. The Dl geometry resulting from such an orbital scan is a thin ring 
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around the patient. Unfortunately, a thin ring geometry for Dl does not produce 
sufficient data for the FBP algorithms in an embodiment of the invention. In an 
embodiment of the invention, a simple linear scan down the axis of the patient body does 
provide sufficient information. This counterintuitive result requires careful explanation. 

[102] The important idea of FBP algorithms is that a 3D delta function in the spatial frequency 
domain permits the identification of the Fourier transform of a backprojected image with 
the Fourier transform of the source distribution. In the preceding derivation two 
functions M and D were evaluated. The function M produced a 2D delta function by 
careful selection of the backprojection weighting; whereas, the function D produced a ID 
delta function from the detector geometry. Without both delta functions, the FBP strategy 
generally fails (spherical geometry is an exception to this rule). For the thin ring 
geometry, no delta function appears in the function D (see Table 1, Geometry 3), so the 
FBP strategy fails. Nonetheless, one might attempt analytic inversion for the thin ring 
geometry. After all, the 2D delta function from the function M yields the ID convolution 
integral of equation (36) that should be invertible by Fourier techniques. This approach 
requires two additional 3D Fourier transforms and an additional filter function. 
Unfortunately, the additional filter generally entails dividing by a (Bessel) function that 
has zeros. The divergences caused by these zeros can only be removed by orbiting the 
camera at two (or three) different radii and taking appropriate combinations of filters. 
One concludes that such a reconstruction may be possible with two or three rings and 
with twice the computer execution time as the FBP algorithms. In contrast, a linear scan 
along the z axis yields a function D with a delta function, so that a FBP algorithm follows 
immediately. 

[103] A Compton camera with simple linear Dl geometry provides similar data. Of course, 
two linear scans along the axis would be better than one, and a cylindrical detector would 
be better yet. However, the preceding analytic analysis demonstrates conclusively that a 
simple ring detector (or equivalently a small primary detector orbiting the patient) may 
not be a good configuration. The lesson is simple: the geometry of the primary Dl 
detector (or equivalently, the surface scanned by the Dl detector) is important in the 
reconstruction. A naiv6 choice of this geometry can doom the best detector electronics 
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and engineering. (On the other hand, the geometry of the secondary detector D2 is 
relatively unimportant. If the D2 detector samples n on a half sphere surrounding each 
point on Dl, then sufficient information will be available for FBP.) 

[104] Figure 9 shows an apparatus 900 for reconstructing a three-dimensional image 
representing a source distribution in accordance with an embodiment of the invention. 
Apparatus 900 comprises Compton camera 100, processor 901, input module 903, output 
module 905, and positioning module 907. Compton camera 100 interfaces with processor 
901 and collects observed data. Processor 901 processes the observed data to reconstruct 
a three-dimensional image of a radiopharmaceutical source distribution in accordance 
with process 600 as shown in Figure 6. 

[105] The reconstructed image may be displayed to a clinician on output module 905. In 
addition, the clinician may provide processing parameters (e.g., index p) to processor 901 
when reconstructing the three-dimensional image. Also, the clinician may instruct 
Compton camera 100 to be repositioned through positioning module 907. Alternatively, 
repositioning module 907 may automatically reposition Compton camera 100 such as by 
linearly scanning along a patient's axis in order to provide a three-dimensional image. 

APPENDIX A: Physical Justification of the Compton Camera Equation (10) 

[106] The integral equation (10) represents a model for Compton camera imaging. In this 
appendix, the physical assumptions underlying that equation are elucidated. The original 
data from a coincident observation in the Compton camera consists of four parameters 
( z , E { ;w,E 2 ) where z6Dl and w6D2. The original data is converted into the 
canonical coordinates in an embodiment of the invention 

(z,E 1 ;w,E 2 )=>(z,n,a f E) 

by the algebraic relations in equation (8): 




,E = 



Ei + E : 



and a = 



( E 2 + E^ - m e Ej 
(E 1 + E 2 )E 2 



(8) 
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[107] Equation (10) can be best understood by examination of the special case of a point source 
of radiation. If one assumes that 

f(x) = 6 3 (x-X 0 ) (A.1) 
and that the original energy of emission is E, then, for p=0, equation (10) gives 



g(2 ' n '°' B)= 2 Jt ]2-x 0 f 6 



n • ( z - x 0 ) 



-a 



(A.2) 



[108] This expression for g requires physical justification. Because g describes the density of 

observed events in the six dimensional space ( z , ii , a , E ), one expects inverse 

dimensions for each of these parameters. Because the vector z is located on a 2D 

surface, one expects cm 2 for that parameter. The vector n occurs on a 2-sphere and, 

therefore, its density is given in steradians (which is technically dimensionless). 
Likewise, the cosine o is dimensionless, whereas, the energy is assumed fixed (integrated 

over a narrow energy window). Thus, the function g has dimensions cm 2 and gives the 

density of scatter events detected in the primary detector Dl at point z that are 
subsequently observed in the secondary detector with parameters ii , a , and E. [In 

equation (A.2) the term | z - x 0 \~ 2 gives units cm 2 , W ( a , E ) gives the units 

steradians \ and finally, the delta function gives the units cosine \] Heuristically, one 
expects that g will be proportional to a product of five terms: 



where <I> = i 



g a <D p, p 2 p 3 p 4 

r y-vzy flux at z in Dl 
originating from a 
^ point source at X 0 



(Probability of an incident y ray ^ 
undergoing Compton scattering 
in the detector Dl 



(A.3) 
(A.4) 

(A.5) 
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p _ I Probability of the y ray j . ^ 

2 " I being scattered in the direction n I * ' 

p _/ Probability of the scattered y ray | . 

3 " ^ propagating from D 1 to D2 J ' 1 7J 

p _ ( Probability of the scattered y ray \ 
ma ^ 4 "{ being detected in D2 J* lA '*' 

[109] Figure 10 shows a detector Dl flux in accordance with an embodiment of the invention. 
As illustrated in Figure 10, the flux O is proportional to the expression 

» cos cp n m 

4jc| z-x 0 J 

where cos <p is the cosine of the angle of the incident radiation with the normal to the Dl 
detector surface. In particular, from Figure 10, one finds that 

[m-(z-x 0 )l 

cosq> = ^ — r-~ — ~i (A.10) 
|z-*o| 

where m is a unit normal to the Dl surface. As illustrated in Figure 1 1, the probability 
(Pj) of an incident gamma-ray interacting by Compton scattering in the Dl detector is 
given by 

where n is the attenuation coefficient associated with gamma-rays of energy E in the Dl 
material, T is the thickness of the detector, and o Compton and a Tota , are the cross sections 

for Compton and all interactions, respectively. In equation (A.l 1), the first expression in 
curly brackets { } is the probability of a Compton interaction compared with all other 
interactions; whereas, the second expression in square brackets [ ] gives the probability of 
the gamma-ray interacting (by any type of interaction) within the detector. The product of 
these two expressions gives the probability of a Compton interaction within detector Dl. 
The probability (P 2 ) of Compton scattering in the direction n is proportional to the 
differential scattering cross section [as given by W in equation (7)], 
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P 2 a W ( a , E ) 6 



|z-x 0 | 



(A. 12) 



where the delta function assures that the vector n lies on the cone with the appropriate 
cosine a. 

[110] The remaining probabilities, P3 and P4, are both assumed constant. In reality, the 
probability P 3 of the scattered photon propagating from a point in detector Dl to a point 
in detector D2 is extremely complicated. Basically, the only phenomenon inhibiting the 
propagation of the scattered photon from detector Dl to detector D2 is a second 
interaction within the primary (Dl) detector. Whether the scattered photon escapes 
detector Dl depends on the thickness of material it must transverse to escape the detector 
and the energy of the scattered photon. Consequently, the probability of escape depends 
on the depth of interaction within the detector crystal and the direction of scatter. This 
probability can best be studied with Monte Carlo simulations. However, for the purposes 
of the embodiment of the invention (namely, development of reconstruction algorithms), 
the probability P3 will be assumed constant. On the other hand, the probability P4 of 

detection of a scattered photon in detector D2 can be accurately estimated as P 4 « 1 . The 

secondary detector D2 is generally designed with thick crystals so that all incident 
photons are absorbed. Thus, the probability of absorption P4 can be assumed constant 
without loss of generality. 



[Ill] Based on the heuristic model of previous paragraphs, the function g satisfies the 
proportionality: 



P 3 P 4 cj Compton 
8 a 2o Total COS * 



1 - exp 



f_j£Ll 

^ cos qp ; 



W(o,E) 
2*1 2-* 0 I 2 



6 



n • ( z - * Q ) 
|*-*o| 



- a 



.(A.13) 



[112] In equation (A.13), the expression in curly brackets { } is the desired result [equation 
(A.2)]. Of the remaining terms, all are constant except those containing cos 9. Offhand, 

this cos q) dependence suggests that equation (A.2) does not accurately model the 
behavior of Compton cameras. However, Compton cameras are generally designed with 
very thin Dl detectors (so that the scattered photons are not immediately reabsorbed in 
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Dl). As a result, the expression in the exponential of equation (A. 13) is generally 
small (|iT«l), so that 



or 



ex P^ coscpJJ cos<p 2cos 2 <p 
cos J 1- exp(- JjJ^ j 



1- 



2 cos <p 



(A.14) 
(A.15) 



[1 13] If cos qp » p.T , one can approximate 



cos qp 



" ex p(-c5?q>)] " ^ T 



(A.16) 



[114] Thus, the cos qp dependence is only important for extremely obtuse angles (cos qp = 0). 
For all other angles, one can write 



. a P 3 P4°Compton^T / W(o,E) fi 
8 " 2 °Total \ 23C| 2-Xo I 2 



P • ( z - % 0 ) 
|*-*o| 



- a 



(A.17) 



[115] The terms containing cosqp in equation (A. 13) are approximately constant because the 
decrease in flux (cos <p) is offset by an increase in the interactions within the Dl detector 

1 1 - exp ( - pT/cos cp ) J due to the increased path length (as shown in Figure 11). If one 

drops the constant coefficients preceding the curly brackets { } in equation (A.17), the 
resulting expression is just that of equation (A.2). 

[116] If the above approximation fails, the terms containing cos qp will reduce the observed flux 

g — especially for obtuse angles (cos cp « 0). For a fixed source position, the angle will 

be smallest for the Dl pixels closest to the source, whereas, the angle will be most obtuse 
for pixels farthest from the source. This observation implies that, if the approximation 
breaks down, the reduction in flux will affect distant pixels more than near pixels. 
Consequently, one would expect that p<0 would provide a better approximation than p>0. 
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[117] Other researchers have used alternative expressions for Compton camera response. 
Legitimate arguments can be made for different approximations. The relative merits of 
these various formulations are unclear at this time. For purposes of comparison, some of 
these alternative formulations will be expressed in terms of the Compton imaging 
operators (C p ) in an embodiment of the invention. One may denote the response of the 

Compton camera by the function "q" (as referenced in a paper authored by R. Basko, et 
al., "Application of spherical harmonics to image reconstruction for the Compton 
camera," Physics in Medicine and Biology, 43, 887-894, 1998), which (using the notation 
of Figure 1 ) can be expressed as 



q(Pl,P2.8) oc (JconeofldA f(x) 
J Possible 

(A. 18) 



Possible 
Emission 
Points 



[118] One can verify the second proportionality in equation (A. 18) by assuming that the 
function f is given by 

f 0 if the distance between x and PI s R 
f(x) = { (A.19) 
0 if the distance between x and PI > R 

for arbitrary R. In that case, one finds that 



q ( PI , P2 , 9 ) a jtR 2 f 0 sin 9 = xR 2 f 0 ^ l-a 2 , (A.20) 

f RP+l 

whereas, C„[f] « (p ° +l)R p (A.21) 

so that p=l yields the appropriate weighting for the function q. One may use the same 
mathematical expression for the function "X" (as referenced in a paper authored by Cree, 
et al., "Towards Direct Reconstruction from a Gamma Camera Based on Compton 
Scattering," IEEE Transactions on Medical Imaging, MI-13, 398-407, 1994) and, 

therefore, are also dealing with the p=l problem. In both cases the sine term ( J l-o 2 ) 
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has no significant affect on the analysis because it can be absorbed in the definition of g 
in the same manner as W in equation (11). Basically, this term arises because they 
parameterize the scatter angle in terms of G rather than o=cosG [ da = -sin0 d0 ]. 

[119] One may introduce an imaging kernel "k" (as referenced in a paper authored by Tomitani, 
et aL, "Image reconstruction from limited angle Compton camera data," Phys. Med, Biol., 
47, 2129-2145, 2002) that can be written in the notation of this paper as 



which implies that p=2. 

APPENDIX B: Evaluation of the Function M for General Indices (q,p) 

[120] According to the notation of equation (26), the regularized form of M given in equation 
(25) can be written in spherical coordinates as 



k(x,z;9) 




(A.22) 




(B.l) 




[121] The | Q x - | term can found from the generating function of the Legendre polynomials 



[60] as 




(B.2) 



[122] Furthermore, from the addition theorem of spherical harmonics, one finds that 




(B.3) 
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so that 

i r* - „t„(2^T) -t- y " (s » ) t^wr (B") 

[123] The oscillatory Fourier terms in equation (B. 1) can be written as 

* k ir 
exp(2 J rixQ x .v)= J o -^=== J k+W2 ( 2* |v| x ) fc J_ k Y k h ( ^ ) [ Y k h ( fi x ) f (B.5) 

and 

expt-lmsCVco) = J o -^=L J u+1(2 ( 2* 1 3 1 s ) J_ u Y u v ( Q, ) [ Y u v ( )]" .(B.6) 

[124] The integrals over Q x and Q s in equation (B.l) follow immediately from the 
orthogonality of the spherical harmonics 

jjd 2 Q Y™(e,cp)[Y u v (e,cp)]* = 6 nu 6 mv (B.7) 
so that 

tf r r exp 2jti x v • Q T - 2m s S • Q J 

//^//"^ - r^i — A = 

'- 7#7f#r»^(^) W^^W^M") (B.8) 



Y n m (S v )[Y„ m (Q <0 )f 
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[125) By substitution of equation (B.8) into equation (B. 1), one finds that 



M q,p(*.«) = 



&t 2 



V |v||5| R^ + P "=o (2n + 1 j 



.-a J 




dx xP-" 2 J n+l(2 (2rt|v|x) exp 



Wf 231 !" 5 ! 8 ) ex P 



- e 



(B.9) 



- E 



[126] The remaining integrals over x and s differ only in the index of the power (p or q) and, 
therefore, only the integral over x will be evaluated. If ( -n -1 < p < 1 ), no regularization 

c 

is required; one can set e = 0 and find [61] that 



CO 

/ 



dx XP-* J n+U2 (2*|v|x) = 



^|*r 2 r(f-f + i) ' 



(B.10) 



[127] Because the sum in equation (B.9) requires all terms n>0; one concludes that equation 
(B.10) is only useful for ( -1 < p < 1 ). However, the subsequent analysis requires p>l 9 

hence the need for the regularization parameter. The regularization used in this analysis 
was chosen because the integral can be evaluated in terms of well-behaved 
hypergeometric functions. The expression [which is valid for £>0, |v|>0, and 
-1 < p+n] contains four terms involving successively higher powers of e. Writing out 
only the lowest order term, one finds that 



00 

Jdx xP-'" J n+lo (2jt|v|x) expj-e-^L" 



2 3iP+" 2 |v| p+ " 2 



r (2 + 2* 2) F f p-n p+n+1 . j_ j_ 3 . e 4 

r(n_£ +1 ) 23 ( 2 • 2 '4'2-4'- (l67C |,, Rc)2 
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(B.ll) 

[128] Because the limit as e -» 0 of the hypergeometric function is one, the limit of equation 
(B.ll) yields 

, L% Jdx x— v,4*M*H ->[-&] - r(j-f If j ^ 

[129] Thus, the regularization does not effect the preliminary results of equation (B.10); 
however, the range of the index p is extended to p>~l . The application of equation (B. 12) 
in equation (B.9) gives that 

(B.13) 

f 1 1 l2 + 2 + 2j 1 [2 + 2 + 2j A Y m (a)[Y m fQ )T 

n4o( 2n + i) r(f-f + i) r(f-f + i) m ^ nWL nlPtojJ 

for all p>-l and q>-l . 
APPENDIX C: The Selection of Special Cases of p and q 

[130] The astute reader has probably wondered why the indices p and q were introduced and 
their role in the reconstruction algorithms. The technical answer is that the proper choice 
of p and q reduces the function M q p in equation (27) to a delta function that permits easy 

inversion of the integral equation (24). The difficulty in equation (27) is the summation 
over spherical harmonics 

M (9 a) « 2 1 r (f + 2 + 2) r (l + f + i) * ym (o U Y m (3 )T 

M ^'<°) n40(2n + l) r (j_| +1 ) r ( ? _a +1 )">^ Yn ^H Y n^cojj • 

(CI) 

[131) Virtually every paper concerning backprojection of Compton camera data [30-35] 
contains a similar summation. However, because other researchers chose values of p and 
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q at the onset of their calculation, there is no way to adjust the coefficients in the 
summation. The ability to make such adjustments is important because the identity 

JomJ-n^t^v) [Y-(aL)f = 6 2 (fi v -^J (C2) 
simplifies the inversion of integral equation (24). The choice of p and q is, therefore, 
determined by a simple condition: the coefficient 

in equation (CI) must be a non- vanishing constant (that is, independent of n). The 
symmetry of equation (C3) permits the interchange of p and q [ p q ) . 

[1321 For large n, Stirling's formula for the gamma functions yields 

n lim K Piq (n) oc n P^- 2 . (C4) 
From equation (C4), one concludes that p + q £ 2 (C5) 

where the equality ( p + q = 2 ) must hold if K p q is to be a non-vanishing constant. 
[133] The special case p=l/2, q=3/2 yields 

K ,„)- 1 r (M) r (f+f) . F (M) 1 (C6) 
m - M n) -(^^)r( f + |)r( t 4)-(2n + .) r(f + I)-4 «*> 

which provides the desired behavior. Because of the symmetry under exchange 

( p ** q ), the case p=3/2, q=l/2 provides the same behavior. Thus, both the cases p=l/2, 

q=3/2 and p=3/2, q=l/2 yield the delta function of equation (C2). 

Applicability of Embodiments 

[134] As can be appreciated by one skilled in the art, a computer system with an associated 
computer-readable medium containing instructions for controlling the computer system 
can be utilized to implement the exemplary embodiments that are disclosed herein. The 
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computer system may include at least one computer such as a microprocessor, digital 
signal processor, and associated peripheral electronic circuitry. 

[135] While the invention has been described with respect to specific examples including 
presently preferred modes of carrying out the invention, those skilled in the art will 
appreciate that there are numerous variations and permutations of the above described 
systems and techniques that fall within the spirit and scope of the invention as set forth in 
the appended claims. 
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